home *** CD-ROM | disk | FTP | other *** search
/ Australian Personal Computer 2002 November / CD 1 / APC0211D1.ISO / workshop / prog / files / ActivePerl-5.6.1.633-MSWin32.msi / _984bc1392f22f26d82217ad09d02adcd < prev    next >
Encoding:
Text File  |  2001-09-04  |  10.7 KB  |  399 lines
  1. package Math::BigFloat;
  2.  
  3. use Math::BigInt;
  4.  
  5. use Exporter;  # just for use to be happy
  6. @ISA = (Exporter);
  7. $VERSION = '0.02';
  8.  
  9. use overload
  10. '+'    =>    sub {new Math::BigFloat &fadd},
  11. '-'    =>    sub {new Math::BigFloat
  12.                $_[2]? fsub($_[1],${$_[0]}) : fsub(${$_[0]},$_[1])},
  13. '<=>'    =>    sub {$_[2]? fcmp($_[1],${$_[0]}) : fcmp(${$_[0]},$_[1])},
  14. 'cmp'    =>    sub {$_[2]? ($_[1] cmp ${$_[0]}) : (${$_[0]} cmp $_[1])},
  15. '*'    =>    sub {new Math::BigFloat &fmul},
  16. '/'    =>    sub {new Math::BigFloat
  17.                $_[2]? scalar fdiv($_[1],${$_[0]}) :
  18.              scalar fdiv(${$_[0]},$_[1])},
  19. '%'    =>    sub {new Math::BigFloat
  20.                $_[2]? scalar fmod($_[1],${$_[0]}) :
  21.              scalar fmod(${$_[0]},$_[1])},
  22. 'neg'    =>    sub {new Math::BigFloat &fneg},
  23. 'abs'    =>    sub {new Math::BigFloat &fabs},
  24.  
  25. qw(
  26. ""    stringify
  27. 0+    numify)            # Order of arguments unsignificant
  28. ;
  29.  
  30. sub new {
  31.   my ($class) = shift;
  32.   my ($foo) = fnorm(shift);
  33.   bless \$foo, $class;
  34. }
  35.  
  36. sub numify { 0 + "${$_[0]}" }    # Not needed, additional overhead
  37.                 # comparing to direct compilation based on
  38.                 # stringify
  39. sub stringify {
  40.     my $n = ${$_[0]};
  41.  
  42.     my $minus = ($n =~ s/^([+-])// && $1 eq '-');
  43.     $n =~ s/E//;
  44.  
  45.     $n =~ s/([-+]\d+)$//;
  46.  
  47.     my $e = $1;
  48.     my $ln = length($n);
  49.  
  50.     if ( defined $e )
  51.     {
  52.         if ($e > 0) {
  53.         $n .= "0" x $e . '.';
  54.         } elsif (abs($e) < $ln) {
  55.         substr($n, $ln + $e, 0) = '.';
  56.         } else {
  57.         $n = '.' . ("0" x (abs($e) - $ln)) . $n;
  58.         }
  59.     }
  60.     $n = "-$n" if $minus;
  61.  
  62.     # 1 while $n =~ s/(.*\d)(\d\d\d)/$1,$2/;
  63.  
  64.     return $n;
  65. }
  66.  
  67. $div_scale = 40;
  68.  
  69. # Rounding modes one of 'even', 'odd', '+inf', '-inf', 'zero' or 'trunc'.
  70.  
  71. $rnd_mode = 'even';
  72.  
  73. sub fadd; sub fsub; sub fmul; sub fdiv;
  74. sub fneg; sub fabs; sub fcmp;
  75. sub fround; sub ffround;
  76. sub fnorm; sub fsqrt;
  77.  
  78. # Convert a number to canonical string form.
  79. #   Takes something that looks like a number and converts it to
  80. #   the form /^[+-]\d+E[+-]\d+$/.
  81. sub fnorm { #(string) return fnum_str
  82.     local($_) = @_;
  83.     s/\s+//g;                               # strip white space
  84.     no warnings;    # $4 and $5 below might legitimately be undefined
  85.     if (/^([+-]?)(\d*)(\.(\d*))?([Ee]([+-]?\d+))?$/ && "$2$4" ne '') {
  86.     &norm(($1 ? "$1$2$4" : "+$2$4"),(($4 ne '') ? $6-length($4) : $6));
  87.     } else {
  88.     'NaN';
  89.     }
  90. }
  91.  
  92. # normalize number -- for internal use
  93. sub norm { #(mantissa, exponent) return fnum_str
  94.     local($_, $exp) = @_;
  95.     $exp = 0 unless defined $exp;
  96.     if ($_ eq 'NaN') {
  97.     'NaN';
  98.     } else {
  99.     s/^([+-])0+/$1/;                        # strip leading zeros
  100.     if (length($_) == 1) {
  101.         '+0E+0';
  102.     } else {
  103.         $exp += length($1) if (s/(0+)$//);  # strip trailing zeros
  104.         sprintf("%sE%+ld", $_, $exp);
  105.     }
  106.     }
  107. }
  108.  
  109. # negation
  110. sub fneg { #(fnum_str) return fnum_str
  111.     local($_) = fnorm($_[$[]);
  112.     vec($_,0,8) ^= ord('+') ^ ord('-') unless $_ eq '+0E+0'; # flip sign
  113.     s/^H/N/;
  114.     $_;
  115. }
  116.  
  117. # absolute value
  118. sub fabs { #(fnum_str) return fnum_str
  119.     local($_) = fnorm($_[$[]);
  120.     s/^-/+/;                               # mash sign
  121.     $_;
  122. }
  123.  
  124. # multiplication
  125. sub fmul { #(fnum_str, fnum_str) return fnum_str
  126.     local($x,$y) = (fnorm($_[$[]),fnorm($_[$[+1]));
  127.     if ($x eq 'NaN' || $y eq 'NaN') {
  128.     'NaN';
  129.     } else {
  130.     local($xm,$xe) = split('E',$x);
  131.     local($ym,$ye) = split('E',$y);
  132.     &norm(Math::BigInt::bmul($xm,$ym),$xe+$ye);
  133.     }
  134. }
  135.  
  136. # addition
  137. sub fadd { #(fnum_str, fnum_str) return fnum_str
  138.     local($x,$y) = (fnorm($_[$[]),fnorm($_[$[+1]));
  139.     if ($x eq 'NaN' || $y eq 'NaN') {
  140.     'NaN';
  141.     } else {
  142.     local($xm,$xe) = split('E',$x);
  143.     local($ym,$ye) = split('E',$y);
  144.     ($xm,$xe,$ym,$ye) = ($ym,$ye,$xm,$xe) if ($xe < $ye);
  145.     &norm(Math::BigInt::badd($ym,$xm.('0' x ($xe-$ye))),$ye);
  146.     }
  147. }
  148.  
  149. # subtraction
  150. sub fsub { #(fnum_str, fnum_str) return fnum_str
  151.     fadd($_[$[],fneg($_[$[+1]));
  152. }
  153.  
  154. # division
  155. #   args are dividend, divisor, scale (optional)
  156. #   result has at most max(scale, length(dividend), length(divisor)) digits
  157. sub fdiv #(fnum_str, fnum_str[,scale]) return fnum_str
  158. {
  159.     local($x,$y,$scale) = (fnorm($_[$[]),fnorm($_[$[+1]),$_[$[+2]);
  160.     if ($x eq 'NaN' || $y eq 'NaN' || $y eq '+0E+0') {
  161.     'NaN';
  162.     } else {
  163.     local($xm,$xe) = split('E',$x);
  164.     local($ym,$ye) = split('E',$y);
  165.     $scale = $div_scale if (!$scale);
  166.     $scale = length($xm)-1 if (length($xm)-1 > $scale);
  167.     $scale = length($ym)-1 if (length($ym)-1 > $scale);
  168.     $scale = $scale + length($ym) - length($xm);
  169.     &norm(&round(Math::BigInt::bdiv($xm.('0' x $scale),$ym),
  170.             Math::BigInt::babs($ym)),
  171.         $xe-$ye-$scale);
  172.     }
  173. }
  174.  
  175. # modular division
  176. #   args are dividend, divisor
  177. sub fmod #(fnum_str, fnum_str) return fnum_str
  178. {
  179.     local($x,$y) = (fnorm($_[$[]),fnorm($_[$[+1]));
  180.     if ($x eq 'NaN' || $y eq 'NaN' || $y eq '+0E+0') {
  181.     'NaN';
  182.     } else {
  183.     local($xm,$xe) = split('E',$x);
  184.     local($ym,$ye) = split('E',$y);
  185.     if ( $xe < $ye )
  186.     {
  187.         $ym .= ('0' x ($ye-$xe));
  188.     }
  189.     else
  190.     {
  191.         $xm .= ('0' x ($xe-$ye));
  192.     }
  193.     &norm(Math::BigInt::bmod($xm,$ym));
  194.     }
  195. }
  196. # round int $q based on fraction $r/$base using $rnd_mode
  197. sub round { #(int_str, int_str, int_str) return int_str
  198.     local($q,$r,$base) = @_;
  199.     if ($q eq 'NaN' || $r eq 'NaN') {
  200.     'NaN';
  201.     } elsif ($rnd_mode eq 'trunc') {
  202.     $q;                         # just truncate
  203.     } else {
  204.     local($cmp) = Math::BigInt::bcmp(Math::BigInt::bmul($r,'+2'),$base);
  205.     if ( $cmp < 0 ||
  206.          ($cmp == 0 &&                                          (
  207.            ($rnd_mode eq 'zero'                            ) ||
  208.            ($rnd_mode eq '-inf' && (substr($q,$[,1) eq '+')) ||
  209.            ($rnd_mode eq '+inf' && (substr($q,$[,1) eq '-')) ||
  210.            ($rnd_mode eq 'even' && $q =~ /[24680]$/        ) ||
  211.            ($rnd_mode eq 'odd'  && $q =~ /[13579]$/        )    )
  212.           )
  213.         ) {
  214.         $q;                     # round down
  215.     } else {
  216.         Math::BigInt::badd($q, ((substr($q,$[,1) eq '-') ? '-1' : '+1'));
  217.                     # round up
  218.     }
  219.     }
  220. }
  221.  
  222. # round the mantissa of $x to $scale digits
  223. sub fround { #(fnum_str, scale) return fnum_str
  224.     local($x,$scale) = (fnorm($_[$[]),$_[$[+1]);
  225.     if ($x eq 'NaN' || $scale <= 0) {
  226.     $x;
  227.     } else {
  228.     local($xm,$xe) = split('E',$x);
  229.     if (length($xm)-1 <= $scale) {
  230.         $x;
  231.     } else {
  232.         &norm(&round(substr($xm,$[,$scale+1),
  233.              "+0".substr($xm,$[+$scale+1),"+1"."0" x length(substr($xm,$[+$scale+1))),
  234.           $xe+length($xm)-$scale-1);
  235.     }
  236.     }
  237. }
  238.  
  239. # round $x at the 10 to the $scale digit place
  240. sub ffround { #(fnum_str, scale) return fnum_str
  241.     local($x,$scale) = (fnorm($_[$[]),$_[$[+1]);
  242.     if ($x eq 'NaN') {
  243.     'NaN';
  244.     } else {
  245.     local($xm,$xe) = split('E',$x);
  246.     if ($xe >= $scale) {
  247.         $x;
  248.     } else {
  249.         $xe = length($xm)+$xe-$scale;
  250.         if ($xe < 1) {
  251.         '+0E+0';
  252.         } elsif ($xe == 1) {
  253.         # The first substr preserves the sign, passing a non-
  254.         # normalized "-0" to &round when rounding -0.006 (for
  255.         # example), purely so &round won't lose the sign.
  256.         &norm(&round(substr($xm,$[,1).'0',
  257.               "+0".substr($xm,$[+1),
  258.               "+1"."0" x length(substr($xm,$[+1))), $scale);
  259.         } else {
  260.         &norm(&round(substr($xm,$[,$xe),
  261.               "+0".substr($xm,$[+$xe),
  262.               "+1"."0" x length(substr($xm,$[+$xe))), $scale);
  263.         }
  264.     }
  265.     }
  266. }
  267.  
  268. # compare 2 values returns one of undef, <0, =0, >0
  269. #   returns undef if either or both input value are not numbers
  270. sub fcmp #(fnum_str, fnum_str) return cond_code
  271. {
  272.     local($x, $y) = (fnorm($_[$[]),fnorm($_[$[+1]));
  273.     if ($x eq "NaN" || $y eq "NaN") {
  274.     undef;
  275.     } else {
  276.     local($xm,$xe,$ym,$ye) = split('E', $x."E$y");
  277.     if ($xm eq '+0' || $ym eq '+0') {
  278.         return $xm <=> $ym;
  279.     }
  280.     if ( $xe < $ye )    # adjust the exponents to be equal
  281.     {
  282.         $ym .= '0' x ($ye - $xe);
  283.         $ye = $xe;
  284.     }
  285.     elsif ( $ye < $xe )    # same here
  286.     {
  287.         $xm .= '0' x ($xe - $ye);
  288.         $xe = $ye;
  289.     }
  290.     return Math::BigInt::cmp($xm,$ym);
  291.     }
  292. }
  293.  
  294. # square root by Newtons method.
  295. sub fsqrt { #(fnum_str[, scale]) return fnum_str
  296.     local($x, $scale) = (fnorm($_[$[]), $_[$[+1]);
  297.     if ($x eq 'NaN' || $x =~ /^-/) {
  298.     'NaN';
  299.     } elsif ($x eq '+0E+0') {
  300.     '+0E+0';
  301.     } else {
  302.     local($xm, $xe) = split('E',$x);
  303.     $scale = $div_scale if (!$scale);
  304.     $scale = length($xm)-1 if ($scale < length($xm)-1);
  305.     local($gs, $guess) = (1, sprintf("1E%+d", (length($xm)+$xe-1)/2));
  306.     while ($gs < 2*$scale) {
  307.         $guess = fmul(fadd($guess,fdiv($x,$guess,$gs*2)),".5");
  308.         $gs *= 2;
  309.     }
  310.     new Math::BigFloat &fround($guess, $scale);
  311.     }
  312. }
  313.  
  314. 1;
  315. __END__
  316.  
  317. =head1 NAME
  318.  
  319. Math::BigFloat - Arbitrary length float math package
  320.  
  321. =head1 SYNOPSIS
  322.  
  323.   use Math::BigFloat;
  324.   $f = Math::BigFloat->new($string);
  325.  
  326.   $f->fadd(NSTR) return NSTR            addition
  327.   $f->fsub(NSTR) return NSTR            subtraction
  328.   $f->fmul(NSTR) return NSTR            multiplication
  329.   $f->fdiv(NSTR[,SCALE]) returns NSTR   division to SCALE places
  330.   $f->fmod(NSTR) returns NSTR           modular remainder
  331.   $f->fneg() return NSTR                negation
  332.   $f->fabs() return NSTR                absolute value
  333.   $f->fcmp(NSTR) return CODE            compare undef,<0,=0,>0
  334.   $f->fround(SCALE) return NSTR         round to SCALE digits
  335.   $f->ffround(SCALE) return NSTR        round at SCALEth place
  336.   $f->fnorm() return (NSTR)             normalize
  337.   $f->fsqrt([SCALE]) return NSTR        sqrt to SCALE places
  338.  
  339. =head1 DESCRIPTION
  340.  
  341. All basic math operations are overloaded if you declare your big
  342. floats as
  343.  
  344.     $float = new Math::BigFloat "2.123123123123123123123123123123123";
  345.  
  346. =over 2
  347.  
  348. =item number format
  349.  
  350. canonical strings have the form /[+-]\d+E[+-]\d+/ .  Input values can
  351. have embedded whitespace.
  352.  
  353. =item Error returns 'NaN'
  354.  
  355. An input parameter was "Not a Number" or divide by zero or sqrt of
  356. negative number.
  357.  
  358. =item Division is computed to
  359.  
  360. C<max($Math::BigFloat::div_scale,length(dividend)+length(divisor))>
  361. digits by default.
  362. Also used for default sqrt scale.
  363.  
  364. =item Rounding is performed
  365.  
  366. according to the value of
  367. C<$Math::BigFloat::rnd_mode>:
  368.  
  369.   trunc     truncate the value
  370.   zero      round towards 0
  371.   +inf      round towards +infinity (round up)
  372.   -inf      round towards -infinity (round down)
  373.   even      round to the nearest, .5 to the even digit
  374.   odd       round to the nearest, .5 to the odd digit
  375.  
  376. The default is C<even> rounding.
  377.  
  378. =back
  379.  
  380. =head1 BUGS
  381.  
  382. The current version of this module is a preliminary version of the
  383. real thing that is currently (as of perl5.002) under development.
  384.  
  385. The printf subroutine does not use the value of
  386. C<$Math::BigFloat::rnd_mode> when rounding values for printing.
  387. Consequently, the way to print rounded values is
  388. to specify the number of digits both as an
  389. argument to C<ffround> and in the C<%f> printf string,
  390. as follows:
  391.  
  392.   printf "%.3f\n", $bigfloat->ffround(-3);
  393.  
  394. =head1 AUTHOR
  395.  
  396. Mark Biggar
  397. Patches by John Peacock Apr 2001
  398. =cut
  399.